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Queueing network models play an important role during each stage 
of a computer system's life cycle (from initial conception to system 
maturity), where in each stage broadly applicable performance anal- 
ysis tools are needed. This paper presents new results which contrib- 
ute to the foundations of a tool to support performance analysis and 
modeling activities. In dealing with some performance issues, it is 
important to be able to quantify distribution or moment information, 
because these quantities can influence system capacity and service 
and performance measures. It is also important that models include 
the effect of congestion adaptive i/o devices, in a stable and efficient 
manner, for this inclusion can significantly affect the outcome of 
studying certain performance issues. We address the problem of 
direct, recursive computation of moments of the queue size distribu- 
tions at a class of service centers embedded in a mixed network of 
queues. The parameterized class includes state-dependent processing 
rates useful in modeling congestion adaptive i/o devices. We also 
present results for calculating moments of both the waiting time and 
virtual delay (work backlog) distributions at a clads of service centers. 
In addition, we obtain a Little's Law type of relation between delay 
moments and queue size factorial moments. For a class of networks, 
an algorithm is given for the direct, recursive computation of the tail 
of the node delay distribution. 

I. INTRODUCTION 

Performance analysis and modeling activities are essential for an- 
swering key questions at various stages of a computer system's life- 
cycle, ranging from initial conception to maintaining and growing a 
mature system. Although both the questions asked and the analysis 
approaches may differ from stage to stage, in each of these stages 
broadly applicable performance analysis tools are needed to support 
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such activities. This paper presents new results which contribute to 
the foundation of one such tool: algorithmic techniques for efficiently 
solving a class of queueing networks. 

In dealing with some performance issues, it is important to be able 
to quantify distribution or moment information (e.g., delay variability 
as opposed to only the mean delay) 1 because these quantities can 
influence system capacity and service and performance measures. It is 
also important that models include the effect of congestion adaptive 
i/o devices, in a stable and efficient manner, for this inclusion can 
significantly affect the outcome of studying certain performance issues 
(e.g., the impact of mmtiprogramming). 2,3 

In this paper, we present results for the direct recursive computation 
of moments of the queue size distributions at a class of service centers 
embedded in a mixed network of queues. The class of service centers 
allows us to efficiently treat, in a stable manner, a parameterized class 
of state-dependent processing rates useful in modeling congestion 
adaptive i/o devices. 3 ' 4,5 By dealing with mixed systems, we allow 
consideration of systems with workloads from a finite population (e.g., 
a collection of terminals), multiprogrammed systems, together with 
workloads from basically infinite customer populations. We present 
results for calculating moments of the waiting time and virtual delay 
distributions at a class of service centers that enable us to quantify the 
variability of node delays, as well as work backlogs. 

The well-known class of multiple resource models, usually referred 
to as product form queueing networks, 6 have been used to address a 
wide range of performance issues, such as capacity estimation and 
planning, bottleneck identification, performance prediction, 7 memory 
interference, and software lockout in multiprocessor systems. 8 Usually 
the models used to address these issues have approximately included 
the factors of interest (e.g., priority processor scheduling disciplines). 9 
A considerable amount of effort has been devoted to the study of this 
class of queueing networks, and efficient computational algorithms 
exist 1,8,10 that allow one to obtain, for example, mean values of the 
desired quantities. While, in principle, the entire network queue size 
state description can be obtained from the above, one may be inter- 
ested in obtaining results, directly, for a more moderate level of detail, 
e.g., moments of queue sizes, as well as in quantifying variability of 
delays at a network node. Existing algorithms for calculating even only 
mean values can become much more complex and sometimes exhibit 
chaotic behavior. 11 This situation arises, for example, when a state- 
dependent service rate is used to model a class of devices, such as an 
efficiently scheduled disk or drum, 3,4,5 * whose efficiency is a function 



* See Ref. 5 for examples of secondary storage units that employ a scheduling 
algorithm which attempts to niinimize rotational latency and/or seek times. 
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of the number of queued requests. The complexity results from a 
requirement that the entire* marginal queue size distribution is needed 
at each step of the computation. 1 

The results in this paper overcome some deficiencies in present 
methods for analyzing networks of queues and enable us to quantify 
important performance measures for which, previously, no efficient 
computational methods existed. In Section 1.1, we discuss how a tool 
for evaluating networks of queues can fit into a performance analysis 
and modeling methodology, and in Section 1.2, we give a more specific 
definition of the problems treated and outline the remaining sections 
of this paper. 

1. 1 Use of network of queues models 

To illustrate the aforementioned need for broadly applicable tools 
and how a queueing network analysis tool can fit into a performance 
analysis and modeling methodology, we consider the various stages of 
a system's life cycle. During the system's conception stage, where one 
is concerned with the services to be offered, broad objectives — per- 
formance, reliability, etc. — basic architecture, proposed components 
and initial sizing, questions involving initial feasibility arise. Initial 
feasibility studies generally require several tools to address the follow- 
ing type of question: Given a proposed architecture and assumptions 
concerning the way the system is to be used (e.g., obtained from a 
gross workload characterization tool which may indicate various usage 
scenarios), how well can the system be expected to perform? The 
answer to this question often leads to a modification of the proposed 
architecture and/or what the system is planned to do. During this 
stage, where system specification is often at a macroscopic level of 
detail, a tool based on a network-of-queues methodology can be useful 
for the performance prediction. Here, for example, use of a tool 
incorporating a central-server 12 queueing network model can be helpful 
in answering questions such as the effects of cpu size, number of disks, 
and the number of terminals that can be supported, while meeting 
broad performance objectives, as well as estimates of cost. The partic- 
ular form of input that may be available from a workload characteri- 
zation is particularly suitable for use by a network-of-queues method- 
ology at this stage. 

During a system design phase (when one is attempting design 
optimization), one may tend to focus more on subsystems, and models 
would tend to include more microscopic details, such as a disk schedule 
model or a CPU process schedule model. Although this stage generally 
requires more detail than is attainable with a network-of-queues model, 



* Some simplifications result when the state dependence disappears above a given 
loading, 110 the so-called limited queue-dependent server. 
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issues such as the effect of the degree of multiprogramming* 2 and the 
effect of a disk schedule which uses a SCAN + algorithm versus a FIFO 
discipline, can be treated. We note that at this stage it is important to 
have models which include, at least parametrically, the effect of a 
device being congestion adaptive since this can significantly influence 
the performance that can be expected.* In this stage, distributional 
information can also be important, as opposed to just dealing with 
mean values, since performance criteria and service objectives may be 
in terms of the tail of a distribution. 

During the stage when one is dealing with and maintaining a mature 
system — adding features, growing, doing capacity estimation and plan- 
ning — a tool or methodology for viewing the overall system is desirable. 
Here, measurement tools, providing resource utilizations, workload 
characterizations, and response times, together with proposed changes 
in system use and growth, can be used with an overall system perform- 
ance model to answer such questions as Is the system adequate? How 
much more load can it handle? or, What system changes are required 
to handle a further load increase? A capacity planning methodology 
which accepts, say, maximum allowable utilizations on various system 
resources and measurements of current utilizations and workloads, and 
predicts allowable increases in load could use the output of a network- 
of-queues tool to determine what the maximum resource utilizations 
should be. We note that the determination of these levels, and thus 
the system capacity, could very well be based on distribution infor- 
mation. 

1.2 Outline of paper and summary 

In this paper we address the problem of recursively computing 
moments of the queue size distribution at a service center embedded 
in a class of product form networks. Parameterizing the service center 
with respect to its state dependence allows us to treat, in an efficient 
manner, congestion adaptive models with either improved or degraded § 
efficiency. We treat mixed systems where the population corresponding 
to a given job type may be fixed (i.e., a closed chain) and where system 
requests corresponding to other job types may arrive exogenously from 
basically an infinite population (i.e., an open chain). Closed chains 



* Reference 2 also shows that the multiprogramming effect, representing an increase 
in potential throughput of up to 75 percent, can depend strongly on the disk scheduling 
algorithm. 

f Actually Ref. 2 considers the LOOK algorithm which is similar to SCAN. 

* Reference 3, which presents a performance comparison of two i/o access disciplines, 
notes that the relative comparisons, representing capacity improvements in excess of 
100 percent, can depend strongly on whether the i/o device is congestion adaptive (e.g., 
a fixed head drum employing the shortest access time first schedule) or congestion 
independent (e.g., fifo). 

8 Resulting, for example, from increased overhead with increased processor queueing. 
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arise when dealing with finite population models where, for example, 
a finite number of terminals places requests into a computer system. 
They also arise in modeling multiprogrammed systems where the size 
of the finite population corresponds to the degree of multiprogram- 
ming. 

Figure 1 shows an example of a mixed system where the closed 
chains correspond to each of the terminal groups, each group with 
possibly different think-time distributions and routing and service 
requirements (workloads). Figure 2 shows a closed network model 
representing requests arriving to a system, S, over a finite collection of 
access trunk groups. A request is blocked if all trunks in the group are 
occupied. Note that a trunk is held during the entire time the request 
is in S. When a trunk is available in group i, requests enter S from 
group i at rate A, (the Poisson arrival rate to group i). The population 
of chain i is Kt, the size of trunk group i, and the blocking probability 
at group i corresponds to 1— utilization of the node with service rate 
A,. In addition to throughputs and resource utilizations as the usual 
quantities of interest, other quantities of interest may be queue size 
moments, and the distribution, or moments, of the time from trunk 
activation to its first entry into the CPU (reaction time). 

We treat the problem of recursively calculating moments of the 
waiting time and virtual delay distributions at first-come-first-served 
(fcfs) state-independent service centers, which enables us to get a 
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handle on the variability of nodal delays and to estimate the distribu- 
tion of such delays. With knowledge of the virtual delays, we can 
compare performance, for example, as viewed by arriving customers to 
that viewed by an outside observer. We can thus estimate performance 
measures such as the system reaction time 13 distribution, which is a 
measure of the time between a terminal request and this request first 
getting the attention of the CPU. In another application, it could 
correspond to the time between a request coming into a system over 
an access trunk group and the start of processing. 

In addition to obtaining algorithmic results, we obtain a relation 
between moments of node delays and queue sizes at a class of fcfs 
service centers. We finish with an algorithm for the recursive compu- 
tation of the tail of the nodal delay distribution at a fcfs service center 
embedded in a network consisting of either single-server, state-inde- 
pendent nodes or infinite server nodes, e.g., the central server model 
of multiprogramming. 12 

The organization of this paper is as follows: In Section II, we present 
the class of networks under consideration, including specification of 
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the job-type characteristics and specification of the class of service 
centers. In Section III, we consider closed systems and present results 
for the factorial moments of queue size distributions, for different 
levels of customer aggregations. The aggregations considered are the 
total number of jobs of a given type at a network node and total 
number of jobs of all types at a network node, the former requiring 
consideration of joint factorial moments and correlations. Results for 
mixed systems are given in Section IV by considering a mixed system 
as the limiting case of a multi-job-type closed system, formed by 
augmentation, as the population increases and the augmented node 
becomes the bottleneck. 

The recursions for nodal delay and flow time moments are presented 
in Section V for closed systems and in Section VI for mixed systems. 
Delay results are obtained as experienced by arriving customers or as 
experienced by an outside observer i.e., the virtual delay. Delay distri- 
bution results appear in Section VII, and the appendices contain 
details of the investigations. 



II. CLASS OF NETWORKS: SERVICE CENTER AND 
JOB CHARACTERISTICS 

We define the structure of the networks in terms of the types of 
service centers or nodes and the job-type characteristics in terms of 
their routing through the network, their service requirements per 
service center visit,* their populations for closed job types, and their 
exogenous arrival rates for open job types. 

We consider a network with R types of customers (jobs, chains) 
where s = 1, 2, • • • , r correspond to closed chains and s — r + 1, 
• • • , R correspond to open chains. The nodes in the network are either 
single-server nodes, using a fcfs, + processor sharing or last-come-first- 
served preemptive resume (lcfs-pr) queueing discipline or infinite 
server nodes. At the single-server nodes, the processing rate can 
depend on the total number of customers present. 

Customers of a given job type (chain) may change class membership 
as they traverse the network but always remain in the same chain. 
Allowing customer class change allows one to model a broad class of 
routing scenarios including a deterministic, fixed sequence of node 
visits and also allows for different visits to, for example, a processor 



* Job characteristics can also be defined in terms of the workload requirements, the 
average resource usage/job lifetime for each resource where job lifetime is suitably 
defined for closed customer types. 

+ Other types of customer selection rules, not depending on actual customer service 
requirements, (e.g., random selection, lcfs nonpreemptive) can also result in product 
form.' 4 
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sharing node by a given customer to have different service requirement 
distributions. 

For the open chains, s = r + 1, • • • , R, we let A,.o.c, s = rate of 
exogenous arrivals of chain s customers, with class membership c, to 
node i. Exogenous arrivals of customers, of a given chain-class pair, to 
a particular node are assumed to be Poisson.* 

If we let p%uc = the probability that a chain s job completing service 
at node k as a class I customer will next enter node i with class 
membership c, then the actual node-class-chain flow rates A,. c , s satisfy 
the so-called traffic equation, 

Ai,c,s = Aj,o,c,s + Zi Ak,l,sPkl.ic', 

s = r + 1, . . . , R; (i, c) e J s , (1) 

where I s is the set of feasible node-class pairs for chains s jobs. 

We note that for each closed chain, the node-class-chain rates satisfy 

Ai,c,s = 2j ^k,l,sP%l,ic', 
(*,/)£/. 

s=l, ...,r; (i, c)Gl 8 . (2) 

The closed chains are further specified by their population, 

K=(#i,tf 2 , -■>,&), (3) 

where K s denotes the system population corresponding to the closed 
chain s. We denote the rate of customers of a given chain, s, flowing 
into a node as 

Ai, s = 2j Ai,c,s> (4) 

rtECi(s) 

where 

G(s) = [c;(i, c)G/ s ]. 

We assume that for each open chain the traffic equation (1) has a 
unique solution + for the traffic rates 

(\ ics ; (i, c) E Is); s = r+l, ••• ,R, 

and that every closed chain has an irreducible routing matrix so that 



* This does not preclude having a large class of state-dependent arrival processes, 
including examples where customer sources are turned off or blocked when their system 
population reaches a threshold. These can easily be transformed into a mixed network 
of the type being considered. 

* This precludes certain pathological cases where customers enter the system and 



never exit. 
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the solution to (2) is unique up to a scalar.* Thus, the actual flow rates, 
for a given job type, can be specified by an arbitrary solution to (2), 
yielding relative arrival rates for classes in a closed chain, and a 
proportionality constant. 

We note that the actual proportionality constant, for a given closed 
chain, say A S (K), depends on the population vector K and that the 
actual flow rates, X ICS (K), can be written as 

X ICS (K) = X S (K)A ICS , (i, c) £ /.. (5) 

Thus, the actual node-chain flow rates Xjs(K) can be written as 

X 1S (K) = XsflOA*, (6) 

where A JS is given by (4). 

The assumptions on the service time distributions are those for 
which the product-form solution holds (see Ref 6). We denote 

IMcl = mean amount of service required by a chain s class c customer 
at node i when only one customer is present at node i (note that 
at a fcfs node it is required that ju, cs = /!,•«>/• for all c, s, e and /), 

and 

Hi(k) = the processing or service rate of node i when a total of k 
customers (regardless of type) are present. + 

Clearly, the average service requirement of an arbitrary, with respect 
to class, chain s customer, denoted nZ\ is given by the weighted 
average of the individual class average service requirements. 

M* 1 - S ^/^. (7)* 

ceC(s) his 

While we allow the network to contain nodes with general state- 
dependent processing rates, we focus on getting queue size moments 
at nodes with a parameterized class of processing rates defined by 



* This follows from the assumption of irreducibility of the routing matrix 

P" = [Pk.c]\ kl, ic E L 

and precludes the pathological case where the set of communicating node-class pairs 
can be decomposed into disjoint subsets. We note that the matrix can be made 
irreducible by using enough chains. 

f For example, the number of seconds of processing per second of elapsed time. 
Without loss of generally we have let /*,(1) = 1. 

* Recall that these are not adjusted for state-dependent processing rates. Another 
interpretation of (7) is that \iu represents the average service time of an arrival to node 
i in a system whose only customer is a single chain s customer. 
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Fig. 3 — A class of state-dependent service rates. 



Mk) = 



CliO + di\k 



(8)* 



where an > 0; a* + an = L We note that Oa = 1 corresponds to a 
single server with state-independent processing rate, while an = 
corresponds to an infinite-server node, often useful in modeling a finite 
population source— a collection of terminals— or a random delay. 
These are shown in Fig. 3, along with cases < an < 1 and an > 1 that 
can be used to represent congestion adaptive devices. We also note 
that the family of service rates may also be useful in approximating 
subsystems with restricted entry 1516 by a single state-dependent node 
as is done in mmtiprogramming models. 17 What is needed here is the 
subsystem throughput as a function of the sizes of its populations, and 
recursive methods are naturally suitable for obtaining this. 

III. QUEUE SIZE FACTORIAL MOMENTS— CLOSED SYSTEMS {R = r) 

We present recursions for two different levels of customer aggrega- 
tion: ( i) for the total number of jobs of all types at a network node, and 
(ii) for the total number of jobs of a given type, the latter requiring 
consideration of joint factorial moments and correlations. 

We denote 



* The methods of this paper can be extended to include a generalization of this family 
of processing rates. The above family provides a good fit to the empirically obtained (by 
R. J. T. Morris) state dependence used in Ref. 2 for an efficiently scheduled (similar to 
SCAN) moving head disk over the multiprogramming range considered. 
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fly(K) = E[ki(ki - 1) • • . (ki -j+ 1); K] (9) 

as they'th factorial moment of the total number of customers, ki, at 
node i at an arbitrary time in equilibrium, for a system with population 
vector K. Appendix A shows that /?«,(K) satisfies the following recur- 



sions 

A(K) = £ ^^ {oaiMK - e.) 

s—l Mis 

+ [1 + (j - l)aa]^-,(K - e.)} , (10)* 

with initialization fi u (0) = 0; j > and /?,„(K) = 1. The quantity e, is 
a unit vector in direction s. We thus obtain the y'th moment at popula- 
tion K by updating the y'th moments corresponding to systems with 
one less customer for each chain in addition to similarly including the 
effect of the (j — l)st moments. The node -chain throughputs, 
A, S (K), are available via standard mean-value analysis; 1 however, the 
standard analysis requires computations + of marginal probabilities at 
nodes with state-dependent processing rates. For networks, with nodes 
of the type under consideration, we can generalize mean- value analysis, 
in a stable manner, without ever computing marginal probabilities and, 
furthermore, obtain the higher order moments. Wheny = 1, (10) yields 
the generalization 

fitiCS) = S ^^ [ctajftiCK - e.) + 1]. (11) 

s-l (t-ia 

Defining the mean node flow time for a chain s customer as the mean 
time a chain s customer spends at node i per visit* to the node i queue, 
Tjs(K), Appendix A shows that 

TUK) = — [1 + OiifiaiK - e s )] (12) 

[lis 

The mean number of chain s customers at node i 

#,i.«(K) - E(kis\ K) 
satisfies 

fiiis(K) = X ia (K)T ia (K), (13) § 

where the (i, s) throughput satisfies [see (6)] 



* In all recursions, terms with a negative population component are zero, e.g., K — e, 
with K = 0. 

f Those computations can become unstable."■" , 

* If a customer is fed back to the same queue after a service completion, as in the 
central server model of multiprogramming, this corresponds to starting a new flow time. 

s The relation with (12) is via Little's Law for chain s customers at node i. 
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1. Initialize /8,y(0) = 0/> 0, /8 /o (0) = 1. 

2. Loop on K until desired population K*. 

3. Compute node flow time means, 7" /s (K), from (12), ieN(s), s = 1,2, . . . ,/?. 

4. Compute throughput proportionality constants, \(K) from (15), s = 1,2, . . . ,/?. 

5. Compute node-chain throughputs, X /s (K), itN(s). s = 1,2 fl, from (14). 

6. Compute 0/ lff (K) from (13), icNfs), s = 1,2 R and 0,-, (K) from (16). 

7. If do not desire fiM) /> 1 or /ff,y s (K) / > 1,go to 2. 

8. If desire Af(K) V > 1, compute jff,y(K) from (10) at desired nodes for 1 </< J. 
If do not want #y s (K) /' > 1,go to 2. 

9. If desire 0f Jt iK*)\ If K = K*-Je s ,initialize Y/ /0/C/S (K) = /3 /c (K), B- 0,1 J. 

If K = K*-ne s n < 4 compute Y/ / j_n,C,s (K) from (21 > as shown in Fi °- 5- G ° t0 2 " 



Fig. 4— Algorithm for queue size moments— closed systems. 

A*(K)=A«A.(K). (14) 

In terms of the flow times, the chain s proportionality constant is 

X s (K)^ — ; s=l,2,---,R, (15) 

where N(s) is the set of nodes visited by chain s customers. To close 
the algorithmic loop, we use 

#i(K)= 2 j8 fl .(K). (16) 



The algorithm (12) -+ (15) -► (14) -»• (13) ^ (16) (see Fig. 4), with 
initial condition /MO) = and with the relative traffic rates, Xu, as in 
the previous section,* represents a simple, stable, + modification of 
mean value analysis for the desired class of state dependencies. With 
the node-chain throughputs from (14), the aggregate higher order 
factorial moments are obtained, in a numerically stable manner, by 
(10). We note that, unlike the mean values, it is only necessary to 
compute higher order moments at those nodes of interest. Thus, at 
the desired node we have the algorithm (see Fig. 4) (12) -» (15) -> 
(14) -> (13) -*■ (16) -* (10). We now turn our attention to obtaining the 
higher order moments for the lower level of aggregation. 

To obtain the higher order (>1) factorial moments of the number of 
each type of customer at a node 

• We recall the relative traffic levels need not be computed for each population vector 
since they do not depend on K. 

+ Note that no subtractions appear in the computations and recall Oa ^ 0. 

720 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1 982 



PuJK) = E[ku(kis - 1) • • • (ku -7 + 1); K], (17) 

we consider the joint factorial moments 
Y..y./. S (K) = E[k ls (k ls - 1) • • • {ku -7 + D 

•£,(£, -1) ... (fe-/+l);K]. (18) 
Note that 

Y..o./. s (K)=/?,./(K), (19) 

the high level /th factorial moment defined by (9) and computable 
from (10), and 

Y..WK) = /?,-,, s (K) (20) 

the desired low level /'th factorial moment. 
From Appendix A, we have the recursion 

yi,>,/,«(K) = — {anYij-i,i+i,s(K ~ e s ) 

IU, 

+ (1 + 2/0,1)7,. y _i./, s (K - e s ) + [1 + (1 + an)(l - 1) 

+ and - l) 2 ]Y«.y-u-i.s(K - e 8 )}; j > 0. (21) 

For/' = the computation is made via (10) and the identification (19). 
We note that for / = 0, /' = 1 (21) + reduces to (13). Making the 
identification (20), gives 

/?, ys (K) = ^^ [fitj-UK - e s ) + anyu-uM ~ ej]. (22) 

(lis 

We note that if we are interested in the Jth factorial moment of the 
number of chain s customers at node i for a target population K*, i.e., 
fiij a (K*), then (21) and (22) are initialized at population K* — Je s with 

Y«.<u„(K* - Je s ) = Pu(K* - Je s ), 

I = 0, 1, • • • , J, (23) 

the simple high-level aggregation result. We further note that (21) and 
(22) only need be updated along parameter direction e s . This is shown 
schematically in Fig. 5 and the steps in the algorithm in Fig. 4. If, for 
example, one is interested in a mean and variance analysis, this 
specializes to 

(3,2s(K*) = ^^- [fi n ,(K* - e s ) + a, mu ,(K* - e s )], (24) 



Terms with a negative subscript are zero. 
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Fig. 5— Schematic for moments of low-level aggregation— closed systems. 



where the correlation is given by 
XUK* - e.) 



Y«.i,i, 8 (K* -e s ) = 



flu 



[a,i# 2 (K* - 2e 8 ) 

+ (1 + 2a l i)/?,i(K* - 2e 8 ) + 1]. (25) 



IV. QUEUE SIZE FACTORIAL MOMENTS— MIXED SYSTEMS 

To obtain the mixed network results, we consider an augmented 
closed system which contains an external node for each open chain. 
The service rate of an external node corresponds to the system arrival 
rate for the corresponding open chain, i.e., Xo« = ^heNs, \i,oy, s = r + 
1, . • • , R. When customers depart in the original system, they are 
routed to the appropriate external node in the augmented system. A 
departure from a given external node, say, corresponding to the open 
chain s, is routed to node i as a class c customer with probability A, 0cs / 
a 0s . Taking the limit of the closed network results as the populations 
of chains r + 1, • • • , R — * °°, and assuming the external nodes become 
the bottlenecks, 1 we obtain the desired recursion for the high-level 
aggregate factorial moments 



* Note that if an internal node becomes a bottleneck, the original mixed system is 
unstable. 
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/3ij(K) = u+(>-i>a„] h. o i g |y _ l( K) 
1 — a.-ip,- 



+ l n UK) o, (*i/MK - e.) 
s _! ]Ui S (l - a.ip") 



+ [! + (/- l)an]/3,j-AK - e.)} , (26) 



where 



rf- £ - (27) 

s— r+1 Pis 

is the unadjusted utilization [unadjusted for state-dependent service 
rates pi(k) ^ 1] at node i, corresponding to all open chains. We note 
that the stability condition is given by 

P?Oa < 1, (28) 

where we note that plan is the limiting (as A, — > co) utilization due to 
customers belonging to open chains. Equation (26) is initialized by 

#.o(K) = 1 
and 

[U0--W MI)[ m 

l — aapi 
the open network factorial moments. Fory = 1, (26) yields 

#i(K) = n P ° = + 2 *'' (K) = [1 + a,iyS,i(K - e.)], (30) 
(1 - aapi) s =i fJ-M ~ anpi) 

with initial condition 

fli(0)= P? „ ■ (3D 

1 — anpi 

The required node-chain throughputs, X IS (K), in (26) can be obtained 
via a standard type of mean analysis. However, as before, we are 
interested in a generalized analysis that does not involve marginal 
probabilities. These throughputs are obtained via the limiting argu- 
ment which yields the recursions 

fituiK) = ^^ [1 + aafiaCS. - e.)]; s < r (32) 

Pis 

for the closed chains and 

fiaJJS) = — [1 + aafiaWl s > r, (33) 

Pis 

for the open chains. We note that in (32), X IS (K), the closed chain s- 
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node i throughput, is to be solved for, whereas in (33), Xu is simply the 
solution to the traffic equation for open chain s. Summing (33) over all 
chains, we obtain 

ff,(K) - : P? o [1 + aafihOQl (34) 

l — anpi 



where 



it 



#i(K)- £ j8a.(K), (35) 

s— r+l 

ffi(K)- J Ai.0K), (36) 

s-l 

which relates the mean value of the total number of open customers 
at node i to the mean value of the total number of closed customers at 
node i. To obtain the desired recursion for fi iu , s < r we use (34) in 
(32), which results in 

fti,(K) ■ k (K) J 1 + anffi(K - e.)], 8 < r; (37) 
and from Little's Law, the mean chain s flow times 

TUK) =—— i -[I + a,i/fr(K - e.)], ■ < r. (38) 

His{l - a,ip») 

The node chain throughputs are calculated from 

X»(K) = A«X«(K), s < r, (39) 

where, X S (K), the chain s proportionality constant (independent of i) 
is given by 

The algorithm, 

(38) -+ (40) -> (39) — (37) -► (36), 
(see Fig. 6) with initial condition 

J&(0) = 0, 
is identical to the algorithm (see Fig. 4) for closed systems, 

(12) -» (15) -> (14) ->(13) -*<M). 

However, the closed chain service rates have been adjusted to account 
for the presence of the open chain customers, i.e., 

^.-/i«(l - (top?). (41) 
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1. Initialize /?,-, (O) from (31) and fijAO) from (29) for 1 < /'< desired order moment. 
PJ, (0) = 0. 

2. Loop on K until desired population. 

3. Compute node flow time means, T /S (K), s< r 
(i.e., for the closed chains) from (38) for itNIs). 

4. Compute throughput proportionality constants X S (K) from (40) for s < r. 

5. Compute node-chain throughputs, \ S (K), s< r, from (39), itNIs). 

6. Compute /9, ls (K) from (37) for s < r, iiN(s) and 0? (K) from (36). 

7. If desire £? (K) or # /1s (K), for some s > r, use (34) for/3?, (K) and (33) with fa (K) 
given by the sum of 0£ (K) and fift (K). 

8. If do not desire /?.(K) / > 1 or jGL(K) / > 1,go to 2. 

f desire fiyiK) J > 1, compute /?-(K) from (26) for 1 < /< J. 
If do not want &.JK) / > 1,go to 2. 

f desire ft iJs (K*) s < r: if K = K* -^.initialize Y/ oes (K) = £ /6 (K). C = 0,1 ^. 

f K = K*-ne s , n<J Compute ?/ »_„ g r(K) from (42) as shown in Fig. 5. Go to step 2 

11. At the desired population K*, if desire ftjjAK*) for s>r, initialize 

Y /ots' K *) " Piti**) B "0.1. • • ..A Compute y/^-^e^lK') from (45) 
as shown in Fig. 7. 



9. 



10. 



Fig. 6 — Algorithm for queue size moments — mixed systems. 

We note, however, that the adjustment is not directly related to the 
actual open-chain utilization (as in the state-independent case 11 ), but 
rather to the limiting open-chain utilization. 

The above algorithm, together with (33) and (34) describing the 
open chains, represents a simple, stable modification to a standard 
mean value analysis for mixed systems with the desired family of state 
dependencies. 

With the node-chain throughputs from (39), the aggregate higher 
order factorial moments can be recursively computed in a numerically 
stable manner from (26). We note that, unlike the mean values, it is 
only necessary to compute higher order moments at those nodes of 
interest. Thus, at the desired node we have the algorithm (see Fig. 6) 

(38) -+ (40) -> (39) -> (37) -+ (36) -► (26). 

We now turn our attention to obtaining the higher order moments for 
the lower aggregation level. 

To obtain the higher order factorial moments of the number of each 
type of customer at a node, /?, ys (K) defined by (17), we require the joint 
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factorial moments y„/ s (K) defined in (18), noting the relations (19) and 
(20). Treating the mixed system as a limiting closed system, we obtain 

.__. A,s(K) .„ . 

yuj/E) = {a,iY,,y-i,/+i. s (K - e s ) 

+ (1 + 2Za„)Y.,y-i./, s (K - e s ) + [1 + (1 + a«)(l - 1) 

+ and ~ l) 2 ]y«,y-i,/-i. s (K - e fl )}; / > 0; s < r , (42) 

which is the same form as (21). For I - 0, j * 1 (42) reduces to (32). 
For j = 0, use (19) and (26). Making the identification (20), we get 

ft,, 8 (K)=^[/^-u(K-e s ) 

flu 

+ aayu-uJK ~ e»)]i s - r • < 43 ) 

If we are interested in /?,„/, s (K*); s < r, then (42) and (43) are initialized 
at population K* - Je„ with 

Y«.o,/, s (K* - Je s ) = #i(K* - ^e,); Z = 0, 1, ■ ■ • , J, s < r (44) 

the high-level aggregation result. Thus the computation proceeds as in 
the closed system case shown in Fig. 5. 

For the open chains, the limiting system argument results in the 
relation 

yi/t|(K) = — {a»iY<,/-U+i,«(K) 

His 

+ (1 + 2Za l i)Y«.,- Us (K) + [1 + (1 + an){l - 1) 

+ a (1 (Z -l) 2 ]Y..y-i./-i. s (K)} j > 0, s > r. (45) 

At I = 0,y = 1 (45) reduces to (33) and for/ = 0, use (19) and (26). We 
note that (45) is not a recursion in K and can be evaluated as an 
aposteriori computation. Thus, if it is desired to compute /fc„/ |S (K); 
s > r at a fixed desired population vector K*, all that is needed is the 
high-level aggregate /?,/(K*); Z = 0, 1, • • • , J, which starts off the (45) 
computation. This is shown in Fig. 7 and the steps in the algorithm in 
Fig. 6. As in the previous section, it is a straightforward matter to write 
down the expressions needed to do a mean and variance analysis. 

V. DELAY AND FLOW TIME MOMENTS— CLOSED SYSTEMS 

We consider the problem of obtaining recursions for the moments of 
the time a customer spends at a fcfs, single-server, state-independent 
node embedded in a closed network. As earlier, we define the node 
flow time as the length of time an arriving customer, either arriving 
from another node or being fed back from the node output, spends 
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• TO GET fi iJs (K*) = y Uos (K'l ; $ > r 

• DENOTE Yf/ps BY (/,£) 



POPULATION 



I/.8J 



K* 
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K* 



• INITIALIZE WITH HIGH-LEVEL MOMENTS /fa(K*J 

Fig. 7 — Schematic for moments of low-level aggregation mixed systems — open chains. 

until he next exits the server. The node delay is correspondingly that 
portion of the flow time spent in the queue. We also present a Little's 
Law type of relation between moments of the node flow times and 
queue size moments. 

Denoting 7V«(K) as they'th moment of the flow time for an arbitrary 
chain s customer at node i, Appendix B shows that 



7v,(K) = i UK r e -> 

i-i M. 



7V,(K - e.) + i 7Vi. s (K); j > 1, (46) 



where the quantity /i, ' is the common mean service time at node i 
(i.e., jS, = His = pa)* From the initial conditions with respect toy 

T, . S (K) = 1, (47) 

and (46), we have 



Mi 



(48) 



We note that an alternate form of (46), which only involves the flow 
time moments corresponding to the job type of interest, is obtained by 
recognizing that 

T ljl (K-e s ) = T U s(K-e l ). 



* Required at a fcfs node to have a product form solution. Since node i has a state- 
independent service rate, n,(k) = 1. 
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It can similarly be shown (using the approach discussed in Appendix 
B) that the y'th moment of the delay distribution corresponding to a 
chain s customer satisfies 

W ,(K)=E X " <K - e ' ) 



/=] Ml 



and 



• Wtj^K -e.)+ir Wu-xAK); j > 1 , (49) * 



Wl , K) ^ML^^(K-e, + l!^|^ (50) 



fory = 1. 

As with the flow time moments, the delay moments at a given 
population are formed by updating the same order moment for reduced 
populations, together with including the effect of the lower order 
moment at the same population. We, thus, note that y'th moment of 
the delay and node flow time distributions satisfy the same recursions, 
(46) and (49), fory > 1, the difference lying in the boundary equation 
(50), which results in 

WijAei) = 0; y>0. (51) 

It is not necessary, however, to have a separate computation for 
Wi Ja (K) since W iJs is obtainable from T ijs by the relation 

W iJa (K) = 2 ^KK-e.) Ti _ [{K _ es) . j > Q (52)+ 

We note that the distributions of delays, as seen by arriving cus- 
tomers, will usually differ from the work backlog at a node at an 
arbitrary point in time (i.e., the virtual delay) as might be measured 
by an outside observer. If we define 

V u (K)=E(Wi(t);K), (53) 

where Wt(t) is the work present at node i at an arbitrary time in 
equilibrium, as the yth moment of the virtual delay, then it is shown 
in the next section that V«/(K) can be obtained from the moments of 
customer delay using the relations 



* The computation of these delay moments or the preceding flow time moments are 
easily included in the algorithm shown in Fig. 4. 

f Note that fory" = 1, this simply states the known result W,-u(K) = /?n(K - e,)//i,. 
Relations (49), (50), and (52) can be written in a form containing only moments 
corresponding to the job type of interest by using Wi>/(K - e.) = Wi>,(K — ei). 
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ViAK) = I ^^- W ut (K) + 1 Vo-i(K); j > 1 (54)* 
and 

Vh(K) = i ^ (I + W ni (K) V (55) 

Having obtained results both for the queue size factorial moments, 
ftij(K), in Section III and the ordinary moments of the flow time 
distribution 7\, S (K) (see above), we now present a relationship between 
them analogous to a well-known Little's Law type of relationship for 
the M/G/l queue in isolation. Appendix B shows that these quantities 
are related via they'-fold summation 

it p 

ft,(K) = £ ••• E X^KjX^K-es,) 



,-i 



The special, single chain, case R = 1 gives 

&/(#> = ( fi X/(/)W(#-./ + l), (56) 

where A,(Z) is the node i arrival rate when our closed system contains 
/ customers and Ty (K - j + 1) is the y'th moment of the node i flow 
time when our closed system contains K — j + 1 customers. Relation 
(56) is analogous to a known result 19 for an M/G/l queue in isolation 
which states 

ft = XTj. (57) 

Noting that V is replaced by the product of flow rates, each with a 
different population, we see that (56) generalizes (57) in the sense of 
holding for a queue embedded in a closed network with the correspond- 
ing, complicated arrival process. 

VI. DELAY AND FLOW TIME MOMENTS— MIXED SYSTEMS 

To obtain the mixed network results, we consider the augmented 
system as in Section IV. We obtain recursions for the delay and flow 
time moments at a fcfs node with a single state-independent server 



* As pointed out in Section VI, the virtual work results also apply to a limited class 
of processor sharing or lcfs-pr nodes; the limitations being that the service times are 
exponential with fi u = pa (i.e., the same assumptions being made for the fcfs nodes). 
Inclusion of this computation in the closed network algorithm shown in Fig. 4 is 
straightforward. 
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embedded in a mixed network. Results are obtained for delays expe- 
rienced by customers in each of the closed chains and for customers in 
any* of the open chains, thus, yielding moments for the virtual delay 
or unfinished work at a node. 

Denoting J5T 1 as the common mean service time at node i (i.e., p., = 
p» = Paf and p° as the actual utilization at node i due to customers 
belonging to open chains, 

„?- I *?. (58)* 

l-r+\ M» 

we obtain 

r»<K)-i *?. ( ? t ~?f r t .(K-«j 

/_! Hi(l — Pi) 

+ l T,,_ llS (K) + 7P? JVi(K - e fl ); < r, (59) 
p, p,(l-p,) 

where T°,(K) is the jth moment of the flow time experienced by 
customers belonging to any open chain at node i and satisfies 

;=ip,(l- pi) Hi{L — p,) 

The recursion (60) is initialized with 



^-M^^ (0) 'wrm- (61) 

the open network flow time moments, and (59) is initialized with 

Tijs(e s ) = TU0). (62) 

The initialization (62) can be obtained by recognizing that the distri- 
bution of the number of customers found at node i by the single closed 
chain customer is identical to the distribution of the number at node 
i at an arbitrary time point, in equilibrium, in a system with no 
customers belonging to closed chains, i.e., K = 0. Since the distribution 
at an arbitrary point in time is identical with that seen by an arbitrary 
open customer, we have (62). 8 

The results for the moments of the delay distributions can be found 
using a similar argument. For the closed chains, this results in 



• We note that from Ref. 20 it is easy to show that delays experienced by customers 
belonging to different open chains have the same moments. 

f Recall that this condition is required at a fcfs node in order to have a product form 
solution. 

* Since node i has a state-independent processing rate, fu{k) = 1 and p? is the actual 
utilization as indicated. 

'An alternate way of obtaining (62) is by use of (59) at K = e„ use of (61) and 
induction. 
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W/ S (K) = l X " { * *' } W#i(K - e.) + L Wij-xM) 

i=i Hi(i - pi) m 

io", 

M«(l - Pi) 

where W° y (K), the yth moment of the delay experienced by customers 
belonging to the open chains, is also the jth moment of the virtual 
delay or unfinished work at node i at an arbitrary point in time. The 
corresponding result for j = lis. 

Wn.(K) = i ^f "! s . } f i + Wm(K - e.) 
/-i ju,(l - p°) \/x, 

+ -,/' ov ? s<r;y=l. (64) 

The moments of ihe virtual delay, or delays experienced by open 
customers, are given by 

W?j(K)=i *f K) W UI (K) 

i=i m(l - pi ) 



+ -' , Wlj-AK); 7>1, (65) 

Hid — Pi) 



and for j = 1 



W?AK)=i^L-(^ + W n ,(K)) + -f!-l; j-1. (66) 
M M<U _ P« ) \M« / 1 - Pi M« 

We note that for both the open and closed chains, the y'th moment of 
the delay and node flow time distributions satisfy the same recursions, 
(59) and (63), for the closed chains and (60) and (65) for the open 
chains, the difference lying in the boundary equations (64) and (66). 
The initial conditions 

W iu (0) = 0; /<r, 

together with (66) and (65), give 



WW) ==t^—o Wlj-m = j\ 

in 1 - pi 



jMl-p?) 



(67) 



the open network result. Using the same argument as was used for the 
closed flow time moments, we have 

Wij.(e.) = WfAO) (68) 

which can be used to initialize (63) and (64). It is not necessary to have 
a separate computation for Wy a (K) if the flow time moments have 
been computed since they are related by 
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W u .(K)-i UK . e,) 7V,(K-e.) 

l-l Mi 



and 



+ p?T°(K-e s ); s<r (69) 



W|(K) = S ^^ 7V,(K) + P ?T° (K); s > r (70)* 



for the open chains or virtual work. 

To obtain the moments of virtual delay for closed systems, denoted 
by Vy(K), we observe that 

V^(K) = lim W?j(K), 
Pi—° 

which results in the relations given in (66) and (67). 

As a final note here we observe that in addition to applying to the 
fcfs nodes as stated, the virtual work results, W"j, V ijy apply to a 
limited class of processor sharing or lcfs-pr nodes + the limitation 
being that service times are exponential with rates ju^ = flu = & . This 
follows from the insensitivity of the stationary queue size distribution 
to the queueing discipline for this case and the memoryless property 
of the exponential distribution. 

VII. DELAY DISTRIBUTIONS 

We consider a multi-job-type closed* network containing either 
state-independent service centers or infinite server nodes, e.g., central 
server model of multiprogramming, and present results for the tail of 
the delay distribution experienced by a type s job arrival to a fcfs 
service center embedded in a product-form network. Appendix C, 
which uses multidimensional generating functions and the known 
relation between the stationary distributions and those seen by cus- 
tomer arrivals to a node, 21 contains the details of the investigation. 

We denote 

W ls (t;K) = P(d l s>t;K) (71) 

as the probability that a chain s arrival to node i is delayed in excess 
of t in a system with population K. Figure 8 shows a single-chain 
example where the delay shown corresponds to the time between a 



* The inclusion of the results of this section into the mixed network algorithm shown 
in Fig. 6 is straightforward. 

f Or other types not affecting the stationary queue-size distributions, e.g., random 
selection or lcfs nonpreemptive. The identification as actual delay moments is, however, 
lost. 

* It is a straightforward matter to treat the limiting mixed network, as well as a class 
of state-dependent nodes (e.g., multiprocessor nodes) by the techniques presented. 
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I/O PROCESSORS 



K TERMINALS 

<f, -*f-1 SINGLE-SERVER NODES 
NODE M INFINITE-SERVER NODE n u = 
W,(f;Af) "P(d,> f.K) 



MEAN THINK TIME 



Fig. 8 — Closed system illustration — delay distribution. 

terminal request and the time the request first gets the attention of 
the cpu in a system with K terminals. References 22 and 23 contain a 
study of the response time distribution (queueing and service) for a 
single node being fed traffic from a collection of terminals, the classical 
machine repair problem with multiple repairmen. This could corre- 
spond to a multiprocessor version of Fig. 8, but without the i/o 
processors. In Ref. 24 the asymptotic behavior, as the number of 
terminals increases, is studied. We note that the methods described in 
this paper can easily be used to obtain results for the response time 
distribution (delay plus service time) for the above example.* In 
general, the state of the art for obtaining response time distributions 
for multiple resource systems is quite limited. + 

Figure 9 shows a representation of two loosely coupled systems with 
shared mass storage devices modeled as two central server models 
with some shared i/o queues. 29 For i/o requests served on a fcfs basis 
in each i/o queue, the delay distribution of interest shown corresponds 
to the delay in accessing a disk, either dedicated or shared, for each of 
the component systems and to the delay in getting each of the cpu's 
for fcfs scheduling algorithms. We note that the population vector 
here, K = (K\, K2), could correspond to the degree of multiprogram- 
ming for each component system. 

We start our presentation of results with the single-chain case and 
then generalize to the multichain case. For the single-chain case, 
letting 



W i (t;K)=P{d i >t;K), 



(72) 



* The single node being a limited queue-dependent server. 10 

* We note that these response times could involve several visits to a given resource, 
as well as visits to other resources. This is studied in Ref. 25 for the single resource 
queue with feedback. Some other recent work is discussed in Section VIII. 262728 
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Fig. 9 — Two systems sharing mass storage devices. 



it is shown in Appendix C that Wi(t; K) satisfies the Mi order 
recursion 

Wdt; K) = yi(t; K - 1) - £ aj(K - l)Wi(t; K-j), (73) 

7-1 

where Mi < M is the number of single-server nodes in the network. 
The coefficients aj(K - 1) are related to the coefficients of Z J , aj, in 

P{Z) = \[[l--Z )=lajZ J , (74)* 

,=i \ [ii ) j-0 

where A, is the relative arrival rate to node i. 
Clearly, 

ctj=(-l) J l PhPh--- P^ (75) 

l<ii<i2< • • • <ij^M\ 

where p, is the relative utilization Aj/ffc,-. The coefficients of (73) 
clj(K - 1) are obtained from 



* We have arbitrarily labeled nodes 1 through Mi to correspond to the single-server 
nodes. Since these nodes are state-independent service centers, in this section we denote 
the common mean service time at node i as pf 1 , i - 1, 2, • • • , Mi. 
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aj(K-i) = (-iy x h(K-i) 

•p h (K-2) •••p ij {K-j), (76) 
which corresponds to (75) with the actual utilizations 

- PljiK - j) = W£zA (77) . 

evaluated at the appropriate population replacing the relative uti- 
lizations pij. Recall that the actual utilizations are available via 
standard analysis. The forcing function on difference equation (73), 
yt(t\ K — 1), can be obtained recursively (see Appendix C) from 

Mi *> -yrrSiW ( 1+ ^i) fit**- x) (78 > 

and the initial condition 

yt(t;l)-e-*pi(l). (79) 

The quantity a* represents the relative loading on all the infinite- 
server nodes 

M A- 
fl-= 2 -. (80) 

We note that for the network example of Fig. 8, 

yi{t\K)=-g—-p 1 (K)(pr + tLit)yi(t\K-l) t 

where p is the fraction of CPU requests feedback to the terminals, 

mean think time 

r = , 

mean cpu service time 

and int represents the point at which we are evaluating the tail of the 
delay distribution in units of cpu service times. 

For the multi-job-type networks, the tail of the node i delay distri- 
bution as seen by a chain s arrival is shown in Appendix C to satisfy 
the multidimensional recursion 

Wi.it; K) = ydt; K - e.) - I a,,, 2 ....,„(K - e.) 

«K/'i+y'2+- • • +jssMi ) 

■ W is (t; K - y'jei - j 2 e 2 - ... - j R e R ) . (81) 

The quantities aj l ,...j R (K — e s ) are the multidimensional analogues of 
ctjiK - 1) in (76), i.e., they are related to the coefficients of Z\ x • 



* Here \,(K) denotes the actual arrival rate to node i for a system with population K. 
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Z» -■- Z J g,a JX ,... jR ,'w. 



P(Z) 



i-l \ s=l/»is 



(82) 



Analogous to (76), we have 

a,,...,„(K - e.) = (-D ; ' + - "A S {[p., i; i(K - e s ) • ■ • p*, : i(K - e s 

{«p,}eL 

- (/, - l)e a )] • • • • [p iiR .R(K - e. - y'iei 
- • • • JR-ie R -i) • • • p*, ^(K - e s - y'ie, - • 



•■ -O'ie-De*)]}, 
where 

0<SA<Mi (83) 

and {i P q} G L corresponds to 

(in < 121 < ••• < i/,i \ 
I J; i™#i,V, (84) 

llR < l2R< • • • ij R R I 

i.e., cij ) j 2 ,...j R (K — e s ) is a sum of products of utilizations with pi pq -, r 
evaluated at the population (K — e s — jid — j2e 2 — • • • ./r-ier-i — 
(p - l)e r ). The forcing function is obtained recursively from 

(a soo (K) + A, S (K)<) 

— yi(t, K - e s ), 

s^(t),s£^*(iUr = OVr^ *(i), 

(a S00 (K) + X is (K)t)8(K) 

# s S(K-e s ) ■ yv ' 

s E #(i), K r = 0Vr<£ 0t*{i), 6(K - e s ) > + 

0; K r >0 for some r£&*{i), 



ydt;K)=- 



(85) 



with initial condition 

, ( ,e) = ^' e ^'' ; re * (i) (86) 

#(i) denotes the set of chains passing through node i, and &*(i) 
denotes the set of chains which either pass through node i or through 
an infinite-server node. The quantity 

S(K)= 2 , Kr x l (87) 

re* (i) 



( 1+ s)' 



f If 5(K - e.) = and K - e s # [if K - e s = 0, use initial condition (86)] y(t\ K) can 
be computed from the first relation for some s £ ®(i), s E &*(i) corresponding to 
K,>0. 



736 THE BELL SYSTEM TECHNICAL JOURNAL, MAY-JUNE 1 982 



the relative chain r loading on the infinite-server nodes 

M \ 

a m = £ — , (88) 

and the actual chain r infinite-server loading 

d TO (K) = £ — = K(K)a m . (89) 

>=M, + 1 fl/r 

We note that it is only necessary to compute yi(t, K) in the subspace 
spanned by chains passing through node i or any infinite-server node. 
For the special case where M\ = M (e.g., the central server model) 
0t (i) =@*(i) and (85) becomes 

'Hit |K| 



y,(t, K) = 



I, 



where 



p is (K)y,(t,K-e s ); 

As | IV | — 1 

sE@{i), K r = 0Vr£^(i), |K|>1, (90) 
K r > for some r ^ @ {i) , 



|K| = X K r . (91) 



VIII. SUMMARY 

This paper has presented contributions to the foundations of a tool 
to support performance analysis and modeling activities aimed at 
answering some key questions at various stages of a computer system's 
life cycle. The emphasis here has been on presenting easily and 
efficiently computable results for calculating distributional information 
and a stable, efficient method for dealing with congestion adaptive 
devices.^ Mixed systems have been considered to allow us the gener- 
ality of dealing with traffic sources which are fundamentally different 
in their behavior. By obtaining results for different levels of customer 
aggregation, we allow one to consider a macroscopic or more micro- 
scopic level of detail. The virtual delay results allow us to quantify 
differences between service as perceived by an arriving customer and 
that perceived by a measuring device. 

We note that many open questions exist in the areas of obtaining 
results related to the distribution of total time a customer spends in a 
subnetwork consisting of several nodes with feedback (e.g., Fig. 1, 



f In a recent paper, 30 a modification of mean value analysis is introduced to eliminate 
numerical instabilities when dealing with general state-dependent service rates. The 
method involves analysis of complementary systems and evaluation of marginal queue 
size distributions. 
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where the time of interest is the response time).* This is an area of 
active research in the literature (e.g., see Refs. 26 and 28); however, 
the results available are fairly restrictive with respect to the network 
topology* or customer paths, and do not apply to, for example, the 
network of Fig. 1. Approximating the moments of the distribution of 
the overall time to transit a network from the individual node flow 
time moments is one possible approach which would have to be studied 
and evaluated. The problem arises because of statistical dependence 
of a given customer's flow times as he sojourns the network. For open 
Jackson networks, Reiman 33 uses a heavy traffic limit theorem to 
obtain a diffusion approximation for the network sojourn times. An- 
other area of importance relates to the inclusion of priorities in, for 
example, the cpu schedule. We note that an approximation technique 
(based on utilization adjustments) does exist 9 for handling a class of 
priority disciplines and can perform quite satisfactorily in many 
cases.* The approximation is such that it enables us to compute 
performance measures using results in this paper. 

APPENDIX A 

Recursions for Queue Size Factorial Moments 

For closed systems, it is known 1 that the marginal, stationary prob- 
ability distributions satisfy 

<* i m l v X * (K) 

Hi(k) s-l P-is 

.p(ki = k-l;K-e a ); k>0. (92) 

While conceptually, the desired moments could be computed from 
recursively computed marginal probability distributions, it is not rec- 
ommended. Other, more computable, approaches could involve the 
use of generating functions. 35 Our approach is to directly obtain a 
recursive relation for the moments. We use (92), (8), and (9) to obtain 
(10), with the indicated initial conditions in a straightforward manner. 
We obtain the required node-chain throughputs, X, S (K), by considering 
the lower level aggregation 



* Or, for that matter, at a processor sharing node embedded in a general closed 
network. The waiting time distribution for a specific closed network consisting of a single 
processor sharing node fed by a single finite population class is treated in Ref. 31. 

t In Ref. 27 a computational methodology is given for obtaining upper and lower 
bounds where an arbitrary network topology is allowed. When applied to two M/M/l 
queues in tandem (for which the exact solution is known) the upper and lower bounds 
are close; however, over 30,000 states were used in the computation at 80 percent 
occupancy. Computational aspects are presented in Ref. 32. 

* Reference 34, which obtains a closed form solution for a two-node closed network 
with priorities, proposes a criterion under which the approximation technique would be 
expected to perform well. 
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k, = (kn, k i2 , ••• , km), 

where kj, represents the number of chain r customers at node i. 
Denoting A?, as the total number of customers at node i, we can write 
the recursive relation 

p,(k,; K) =-i ** } p,(k, - e s ; K - e s ); ku > 0, (93) 

which is obtainable from the product form solution. From (93), (8), the 
standard Little's law argument at each node and about the entire 
system, and the irreducibility of the routing chains we get (12) through 
(15). To obtain the desired recursion for the joint factorial moments, 
we use (8) and (93) in (18), make the appropriate identification corre- 
sponding to a system with reduced population and obtain (21). 

The mixed system result (26) is obtained directly from (10) by 
considering the augmented system with population 

K' = (K\K r+i , •• ,K R ), 

denoting 

ft(K) = lim fty(K') 

and decomposing the sum over the open and closed chains. The lower 
level aggregation results are obtained by use of 

7 77-T — A(k, - e s ; K); s>r 

fas mKkopLu 

(94) 



lim p,(k,; K') = 



kuHi(ki)ni 



p«(k-ei;K-e.); s<r, 



where e* is an /^-dimensional unit vector in direction s as opposed to 
e s , the corresponding r dimension vector for s < r. 

APPENDIX B 

Recursions for Node Delay and Flow Time Moments 

Denoting W is (t, K) as the complementary delay distribution expe- 
rienced by a chain s customer at node i when the system population is 
K, we can write 

W«(*,K)- I ,7 „. e-WpiM u)>k;K), (95) 

a_i (k — 1)! 

where ju, s = /x,- at the fcfs state-independent node under consideration, 
ki(tT s ) represents the number of customers at node i seen by a chain s 
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arrival to node i, and |K| = Ki + K 2 + • • • + Kr. For the class of 
closed* systems, we are initially considering 

p(ki(tu) = k; K) = p{ki = k; K - e 8 ) , (96) 

i.e., the distribution as seen by an arriving chain s customer is equal to 
the distribution at an arbitrary point in time in equilibrium (i.e., the 
stationary distribution) for a system with one less chain s customer. 20 
Using (96) in (95), we can obtain the Laplace-Stieltjes transform of the 
flow time distribution for node-chain pair (i, s) 

Tu,(v, K) = -£ ^ 'T P(k^k-,K- e.) (-£-) (97) 

The ;th moment, obtained by differentiation, satisfies 



IKI-l 



T iJS (K) = J - Tij-UK) + L I k(k + 1) • • ■ (k + j - 2) 
Hi M« *-i 

• p(ki>k;K-e s ), (98) 

and the summation in (98) can be written as 

l K '- 2 * Xu(K - e 8 ) 
I (g + 1)... (q + j-D 1 

q-o l-l M' 

.p(ki>q;K-ei-e.), (99) 

where we have used (92) to get the one-step recursion on the tail of 
the marginal queue size distribution. Inserting (99) in (98) and making 
the appropriate identification, we get (46), with initial and boundary 
conditions as indicated. 

To obtain the Little's Law type moment relation, we use the defi- 
nition of pij(K) and summation by parts to get 

frj(K)=j e'(/+i)</+2) ••• (l + j-l)p(kt*l + j;K), (100) 

/-o 

which upon repeated application of (92) results in 
iMK) - 2 • • • 1 X^,(K)X JS2 (K - e.,) • • • X l8 , (K - £ e 

.l,* K Z (l+l)...(l + j-l)p(ki>l:K- l e s V (101) 

Hi 1-0 \ 9-1 / 

Identifying the last term, we obtain 



* The mixed system results obtained by the usual limiting system argument are 
reported in Section VI. 
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R R / j-\ 

Pa(K) = I • • • I X JSl (K)A, S2 (K - e.,) ■ • • a*, K - £ e 4 

S, = l S-l \ 7-1 

xTosjIk-I e s \ (102) 

which for the special, single chain, case R = 1 gives the moment 
relation 

MK) = ( II MD J Ty(K - j + 1) . (103) 

V-/f-y + l / 

APPENDIX C 

Recursions for the Tail of the Node Delay Distributions 

We consider a closed system of Mi state-independent, single-server 
nodes and M — Mi infinite-server nodes and obtain an M \ h order 
difference equation for the tail of the customer delay distribution at a 
fcfs node. Denoting 

W is (t; K) = P(dis > t: K) 

as the probability that a chain s arrival to node i is delayed in excess 
of t in a system with population K, we can write 

|K|-1 („./)*-! 

W is (t;K)= X %J—e-"p(ki=»k;K-e.). (104) 

*=i (k — 1)1 

We denote the product-form normalization constant* as 

G(K) = Jg *i(ki)£ 2 (ka) ••• gm(M, 

I 

where 

k/i\ ■ • • kjd \HjiJ \hjr 

>=1,2, ...,M,, 

X \*>' /X * 



»^-)- — ri n — ; y = Mi + i, ...m, (105) 

and we have labeled nodes j = 1, • • • , M\ as the single-server nodes. 
Multiplying (104) by G(K - e s ) and obtaining the generating function, 
we have 



* Our final result will not involve computation of the normalization constant, the 
calculation of which can result in numerical problems. 
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swii^z. 



Li.(Z) = 



R x 

r-l M< 



exp 



■/«,/(i-Z— * 

r-l [H 



, (106) 



where 



and 



LUZ) = I ■ • • £ W*(f; K + e.)G(K)Zf • • • • Zf», (107) 



A",=0 tf„=0 



S,(Z)= S ••■ I ^(K)Zf' ...zJK 

tf,-0 /f H =0 



(108) 



where #,-(K) is the normalization constant for a reduced network 
(node i absent) 

2MK) = £ £i< k ') • ■ ■ ^-i(k l -i)^ + i(k, +1 ) • • • g R (k.R) . (109) 



Inserting 



where 



S,-(Z) = n*(Z), 

Mi 



(110) 



SAD = < 



f 


1 






(l 


R 

-7, 


V 


»y 


\ 




tt* 


/ 


fl 








I 


'^Z r 






gr=lM,r . 







;y<M,, 



j>Mi. 



(in) 



into (106), we get 



/ ff Xir \ I (a r .+X„/)Z, 

P(Z)L ls (Z) = i-* X-Z, r 1 4 Oft Z), (112) 

\r=l M« / 

where P(Z) is the polynomial given by (82) and a m given by (88). 
Inversion of (112) and division by the appropriate normalization con- 
stant results in 



5 *J UK - e s - yid - ... - >fle R ) 

E ■•• 2 «>!.■••>« 

/i-0 yn-o 

Jl+'-'JHSMl 



G(K-e s ) 
hdt,K-e s ) 



Ayi(t,K-e s ), (113) 



G(K-e s ) 
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where «,,,...,„ is the coefficient of Z{ l Z j £ • • • Z j g in the polynomial (82), 
hi(t\ K) is the inverse generating function of Hi(t; Z) given by 



hi(t, K) - e-"'' 



n 

r&dt'ti) 



(a rx + Xirt) 



K r \ 



XirKr 



reSlV) Mi(Oroo + Xirt) 



(114) 



0t(i) denotes the set of chains passing through node i and @*(i) 
denotes the set of chains which either pass through node i or an 
infinite server node. The quantity /, S (K) in (113) satisfies 

UK) = W is (t, K + e s )G(K) . (115) 

As it stands (113) represents an M l \ order recursion for the distribution 
tail; however, the coefficient and forcing function involve the normal- 
ization constant. Upon use of (115) in (113), and recognizing that 

G(K - e s - yiei - ... - j R e R ) 

G(K-e 5 ) 

can be written as a product of node throughput proportionality con- 
stants 



n X 1 (K-e s -/ iei ) 

L*i" 



/,-0 



72-1 

I] X 2 (K - e s - j\ei - I&2) 



I] X R (K- e s -y'iei - 



/ fi =o 



- JR-ieR-i — l R e R ) 



we obtain (81) and (83). 
The recursions (85) for the forcing function 

hi(t, K) 



ydt; K) = 



G(K) 



follow from 



hi(t, K) = 



(a so = + X is t) 5(K) 



K 8 «(K - e . 

where 5(K) is given by (87), and from 

(a soo + \ut) 



hdt, K - e s ); 



(116) 



'<*'), (H7) 



hi(t,K) = 



K 8 



hi(t,K-e s ); s£8(i), sE@*(i). (118) 
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